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ABSTRACT 


A full three-dimensional finite element-multibody structural dynamic solver is coupled to a three- 
dimensional Reynolds-averaged Navier-Stokes solver for the prediction of integrated aeromechanical stresses 
and strains on a rotor blade in forward flight. The objective is to lay the foundations of all major pieces of an 
integrated three-dimensional rotor dynamic analysis — from model construction to aeromechanical solution 
to stress/strain calculation. The primary focus is on the aeromechanical solution. Two types of three- 
dimensional CFD/CSD interfaces are constructed for this purpose with an emphasis on resolving errors 
from geometry mis-match so that initial-stage approximate structural geometries can also be effectively 
analyzed. A three-dimensional structural model is constructed as an approximation to a UH-60A-like fully 
articulated rotor. The aerodynamic model is identical to the UH-60A rotor. For preliminary validation 
measurements from a UH-60A high speed flight is used where CFD coupling is essential to capture the 
advancing side tip transonic effects. The key conclusion is that an integrated aeromechanical analysis is 
indeed possible with three-dimensional structural dynamics but requires a careful description of its geometry 
and discretization of its parts. 


INTRODUCTION 

The objective of this paper is to lay the foundations 
of an integrated three-dimensional (3D) aeromechanical 
analysis for helicopter rotors. The words “integrated 
3D aeromechanics” are used to mean three-dimensional 
structural dynamics coupled to three-dimensional fluid 
dynamics for the prediction of dynamic stresses and 
strains with an equal fidelity of representation in struc- 
tures and fluids. 

The state of the art in aeromechanical analysis uses 
Reynolds Averaged Navier-Stokes (RANS) based Com- 
putational Fluid Dynamics (CFD) containing tens of 
millions of grid points on hundreds of cores, routinely, 
in a research environment for the rotor, and even for 
the entire helicopter. Arbitrary blade contours, air- 
craft configurations, and flight conditions can be an- 
alyzed effectively from first principles. Computations 
in the structural domain continue to use engineering- 
level beam-multibody models that are historically part 
of lifting-line comprehensive codes carried out on a sin- 
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gle processor [1], Thus the high-fidelity CFD airloads 
are not directly integrated with stress/strain calcula- 
tions. Similarly detailed finite element stress analysis is 
performed routinely on structures but only on isolated 
components with loading from beam analysis or previ- 
ous flight test data. Thus high-fidelity structures is not 
integrated with high-fidelity airloads. The need for high- 
fidelity structures for dynamics is clear. Modern rotors 
contain flexible components near the hub that provide 
the critical couplings for dynamics and encounter the 
critical stresses. These have short aspect ratios, open 
sections, and end constraints, and cannot be treated as 
beams. Blades envisioned with advanced shapes and the 
ability to morph in flight cannot be treated from first 
principles using beams. 

The broad objective of this research is to close this 
gap by developing a scalable 3D solid-multibody based 
Computational Structural Dynamics (CSD) solver. The 
formulations for 3D scalability and multibody dynamics 
were covered earlier in Refs [2] and [3]. The focus of 
this paper is on integrated aeromechanics. The specific 
objectives are: 1. to develop a 3D CFD/CSD coupled 
solution procedure and demonstrate it on a realistic ar- 
ticulated rotor and 2. to implement an entire 3D work- 
flow from model construction to aeromechanical solution 
to stress/strain calculation. 
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Figure 1: Integrated 3D workflow. 


The field of computational 3D fluid-structure inter- 
action is vast and varied (see recent reviews and texts 
[4]- [8]) but have never been applied in the context of ro- 
torcraft. Rotorcraft requires a special solution procedure 
— a straight forward exchange of states and accompa- 
nying iterations will not do. Rotor control angles must 
be simultaneously determined to satisfy the aircraft trim 
state. Adequate aerodynamic damping must be provided 
to the structure as a pre-conditioner so that all classes 
of rotors can be solved. Theoretical considerations of 
conservation and consistency are not sufficient, an effec- 
tive analysis must accommodate early-stage structural 
descriptions with significant geometry nris-match. These 
are the particular emphases of this paper. 

The 3D structural model constructed for this work 
is an idealized representation of a fully articulated UH- 
60A-like rotor. Construction of exact geometries from 
CAD and generation of 3D finite element-multibody 
meshes is a separate dedicated effort in itself and de- 
scribed in a companion paper [9]. 

Scope and Organization of Paper 

The main emphasis of this work is on the develop- 
ment and demonstration of an integrated 3D CFD/CSD 
analysis capability. A complete capability demonstra- 
tion also requires all supporting pieces of a new work- 
flow to be covered. Many of these supporting pieces are 
covered in an idealized manner. For example 3D geom- 
etry and meshing is a critical piece for such a capability 
yet dealt in a crude and coarse manner adequate only 
for the purposes of capability demonstration. The com- 
parisons shown at the end between predictions and mea- 
surements are therefore not true validations but only a 
basis for qualitative verification. Nevertheless they pro- 


vide the first glimpse of the kind of detailed stress /strain 
analysis possible using integrated 3D aeromechanics. 

The paper is organized as follows. Following intro- 
duction, the paper begins by describing the new work- 
flow. The next section describes the 3D CSD solver. 
The external CFD solver used for coupling is then briefly 
summarized. The fourth section describes the construc- 
tion of 3D structural models. The fifth section covers 
the 3D CFD/CSD formulation — interfaces (spatial cou- 
pling) and the solution procedure (temporal coupling). 
The final section presents results from a fully integrated 
3D analysis. The paper ends with some concluding ob- 
servations. 

INTEGRATED 3D WORKFLOW 

Integrated 3D analysis requires a new workflow. 
The key pieces of this workflow are envisioned in Fig 1. 

The workflow begins with a CAD model. It can 
be a gross description in the early stages of structural 
design, then progressively refined and populated with 
details as the design advances. The CAD geometry is 
then interpreted into a Structural Analysis Represen- 
tation (SAR) geometry by the dynamicist. The SAR 
includes only those parts considered necessary for anal- 
ysis and identifies them as a flexible part, a multibody 
part, or a device part. Flexible parts are marked for 
3D meshing, multibody parts are assigned their func- 
tionality (joint type and constraints), and device parts 
are noted for special-purpose treatment. The flexible 
parts are those with significant strains. The multibody 
parts are idealizations of constraints that allow arbitrar- 
ily large relative motions between flexible parts. Devices 
are special-purpose parts which allow off-line character- 
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istics to be included (e.g. look-up tables for a non-linear 
lag damper). The next task is to mesh the flexible parts. 
Once meshed, all parts are re-integrated into a Struc- 
tural Analysis Model (SAM) geometry. The SAM task 
merges all 3D parts, prescribes material properties, des- 
ignates constrained nodes for multibody connection, de- 
fines the matematical properties of all joints (rotation or- 
der, locked versus free states and command signals), de- 
fines the physical properties of all joints (stiffness, damp- 
ing, actuation force and ± free play), and completes the 
appropriate description of devices. The completed SAM 
is then ready for the CSD solver. All of these tasks fall 
under the broad category of 3D geometry and meshing 
(a companion paper [9] is devoted to this category). 

A CSD solver must be equipped with special capa- 
bilities for rotorcraft — a pure structural dynamic anal- 
ysis will not do. An internal aerodynamic module is 
needed for efficient coupling with CFD. A trim module is 
needed to achieve the mean rotor operating state. These 
require a top level description of aircraft geometry and 
configuration. Trim requires integrated loads (inertial 
and aerodynamic) at the hub from multiple load paths 
across multiple parts for which a specialized hub loads 
module is needed. A high-fidelity interface is needed for 
3D fluid-structure coupling. The solution procedure for 
coupling is special and requires all of the above modules. 
All of these pieces are grouped within the broad cate- 
gory of a rotor dynamic solver. The 3D fluid-structure 
interface piece is generic to any external CFD solver 
(with structured or unstructured surface mesh). The 
CFD solver used in this work is part of the HPCMP 
CREATE™ - AV Helios software. This required the 
CSD solver to be incorporated within the Helios envi- 
ronment. This environment is presently equipped to ex- 
ecute a CSD solver from only a single processor — as per 
requirements of current generation comprehensive codes. 
While this is not a barrier for the main emphasis of this 
paper it prevents the Structural Partitioning piece of the 
workflow from implementation (Fig 1). It also prevents 
meaningful timing conclusions to be drawn. 

3D ROTOR DYNAMIC SOLVER 

The 3D rotor dynamic solver consists of a core 3D 
solid-multibody CSD solver, an aircraft geometry mod- 
ule, an internal aerodynamics module, a hub loads mod- 
ule, a trim module and a 3D fluid-structure interface 
module. The fluid-structure interface module is de- 
scribed in detail in a later section. The present work 
considers only an isolated rotor so the aircraft geometry 
module is not used. A brief summary of the rest are 
provided below. 

3D Solid-multibody CSD Solver. The flexible parts 
of a structure are solved by 3D governing equations of 
motion derived in a rotating frame (non-rotating is a 
simplification) using generalized Hamilton’s Principle. 


The formulation uses Green-Lagrange strains and second 
Piola-Kirchhoff stresses for strain energy and follows a 
geometrically exact nonlinear Total Lagrangian formu- 
lation. Exact geometry and strains are a pre-requisite 
for unifying multibody dynamics within 3D. The stress- 
strain relationship is linear. Isoparametric, second or- 
der, brick elements: 27-node hexahedral elements and 
10-node tetrahedral elements are available for discretiza- 
tion of the structure. See Ref [2] for details. 

The multibody connections and constraints use Eu- 
ler angle joints. A special formulation is used which 
preserves kinematic exactness between 3D elements un- 
dergoing arbitrary rotations relative each other while at 
the same time eliminating all algebraic constraint equa- 
tions. The special formulation is a joint that is associ- 
ated with 12 states — 6 are interface states that con- 
strain the position and orientation of the joint in space 
and another 6 are the joint states that describe its de- 
formation. The joint states can be constrained (locked), 
commanded (prescribed) or actuated (forced). Connec- 
tions are also special in 3D — they represent true con- 
nections. Connection points must be identified precisely 
on the flexible parts and a mesh generated with nodes 
available at those locations. If the connection points are 
not known a priori a generic connection can be used that 
constrains an entire face of the flexible part. The inter- 
nal stresses at the edge will then be incorrect but the 
kinematic exactness will still be preserved. See Ref [3] 
for details. 

Three classes of algebraic solvers are available: di- 
rect, iterative, and eigen. The direct solver is a skyline 
solver. The iterative solvers are iterative-substructuring 
solvers based on Finite Element Tearing and Intercon- 
necting - Dual Primal preconditioners and equipped with 
Conjugate Gradient and Generalized Minimum Residual 
updates. The eigen solvers are not developed as part of 
this work but available as calls to external scaLapack 
routines [10]. The direct solver is most efficient on a 
single processor but not parallel. The iterative solvers 
are parallel and scalable and meant for large scale dis- 
tributed execution. The eigen solvers are parallel but 
not scalable. Because the present environment restricts 
the solver to a single processor the direct solver is used. 

The algebraic solvers are the building blocks of 
structural solvers. Two types of structural solvers are 
available, both employing a direct time integration of 
the second order nonlinear structural dynamic equations. 
The first is a single time-step Generalized-a Method [11] 
and the other is a two time-step trapezoidal plus three- 
point backward Euler combination [12], The 3D solid- 
multibody models are strongly nonlinear and the dy- 
namic stiffness matrix must be updated at every time 
step for accuracy of the response solution. Additionally, 
structural sub-iterations (Newton-Raphson) are possible 
within both classes of solvers, but required only for rel- 
atively large time steps. 
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Internal Aerodynamics. The internal aerodynamics 
model, presently, is only meant to support CFD cou- 
pling, hence elementary: quasi-steady linear aerodynam- 
ics with uniform inflow. Its purpose is to provide airload 
sensitivities to blade deformations (aerodynamic damp- 
ing) so non-circulatory airloads are also required (for tor- 
sion damping). Aerodynamic angles (angle of attack a 
and sideslip /?) are extracted from the 3D deformation 
field using the deformed chord orientation relative air 
flow and blade velocities at 3/4 — c location (see sec- 
tion on 3D Structural Models) as per thin airfoil the- 
ory. Both angles and rates are required for the angle 
of attack. There is no aerodynamic damping in lag so 
either a damper model is needed (a simple linear model 
is used here) or artificial damping added and bled out 
consistently to decay initial transients. Transients in lag 
dynamics make a direct integration in time to periodic- 
ity an expensive solution process. Even with the known 
damper value for this existing rotor about 15 revolutions 
are needed to decay any lag transients. 

Hub Loads. The root shears and moments in the 
rotating frame are calculated either by force summation 
or from joint reactions. The counterpart of the deflec- 
tion method in beams is a direct stress integration in 3D. 
Stress integration is available for the flexible parts but at 
the ends or at the joint connection points — where the 
integrated loading is actually desired — 3D edge effects 
or local concentration effects occur and contaminate the 
solution. The effects are real but require a high local 
mesh resolution to capture. The contamination of global 
trim solution with local mesh resolution is considered un- 
acceptable. Thus force summation or joint reactions are 
needed instead of the stress integration method. (Note 
that even in beams shear forces cannot be resolved di- 
rectly — if at all — using deformations and a force sum- 
mation is always needed). These are needed also for 
elements that are idealized to be rigid. The force sum- 
mation method is same in principle as in beams except 
that a volume integration is needed in 3D. For a con- 
verged response both should be identical (see section on 
3D CFD/CSD Analysis Results) even though in general 
force summation is more accurate for higher frequen- 
cies and larger time steps. Joint reactions are needed 
for multiple load paths. Here joints are used as sensors 
which means they have to be assigned a small flexibility 
in the direction (and type) of the desired loading. Thus, 
strictly, sensing alters the dynamic characteristics of the 
model — but is not inconsistent with the real behavior 
of systems. 

Trim Solution. Once hub loads are obtained the 
trim solution is straight forward. Presently only isolated 
rotor options are available: wind tunnel trim (targets are 
thrust or coning and two cyclic flapping angles), moment 
trim (targets are thrust or coning and two hub moments) 
and propulsive trim (targets are lift and propulsive force 
and either two cyclic flapping angles or two hub mo- 


ments). The control inputs are via joint commands (ad- 
ditionally shaft tilt for propulsive trim) imposed either 
at the blade root bearing or at the base of the pitch link. 

3D CFD SOLVER 

The CFD solver used is part of the HPCMP 
CREATE™ - AV Helios software. It is an integrated 
capability consisting of an unstructured, node-centered, 
implicit RANS (Spalart-Allmaras turbulence) near-body 
solver, an overset Cartesian, explicit, Euler off-body 
solver, and an implicit hole-cutting based domain con- 
nectivity algorithm. The solver has been extensively val- 
idated with current state of the art beam models for the 
same rotor used in this paper. A description of its gen- 
eral architecture can be found in Ref [13]; validation of 
its aerodynamic and CFD/CSD capabilities can be found 
in Ref [14]. 

The aerodynamic set up is identical to Ref [14]. A 
coarse mesh set-up is used. The near-body grid contains 
4.5 million points and extends up to 5.8 chords. The 
off-body grids contains about 14 million nodes (finest 
spacing of 0.18 chord) and extends up to 58 chords. An 
azimuthal discretization of 0.25° is used with 25 sub- 
iterations for the near-body solver and a single sub-step 
for the off-body solver. The Courant-Friedrichs-Lewy 
number for the off-body solver for this time step is 1.67. 
The CFD solver is executed on 128 processors. 

3D STRUCTURAL MODELS 

A four-bladed fully articulated rotor model is con- 
sidered. Each blade is an assembly of several flexible and 
joint parts. The lag damper is not modeled as a sepa- 
rate device but included as a linear damper as part of a 
joint. Two different configurations (baseline and simple) 
and two different blade meshes (baseline and simple) are 
constructed for this study. The baseline configuration 
and the baseline blade constitute the full-up model. All 
results shown are from the full-up model unless otherwise 
mentioned. 

The baseline configuration consists of four flexible 
parts, three joint parts and two load paths (see Figs 2 
and 3). The flexible parts are the blade, the hub block, 
the pitch horn and the pitch link — all modeled using 
3D elements. The blade contains 592 or 48 hexahedral 
elements (see later), the hub block 8, the pitch horn 3 
and the pitch link 3. Joint 1 is a spherical joint located 
within the hub block and represents the elastomeric flap- 
lag-pitch bearing at the root end. It is connected to three 
flexible parts — the blade, the hub block and the pitch 
horn. It transfers blade loads to the hub block and pitch 
horn control motions to the blade. Joint 2 is another 
spherical joint located between the pitch horn and the 
pitch link. Joint 3 is a slider joint located at the bot- 
tom of the pitch link and represents the connection to 
the swasliplate. It is free to roll and pitch relative to 
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Figure 2: 3D model of an UH-60A -like rotor. 

the swashplate but locked in yaw. The control angles 
are imposed as a linear motion command at this joint. 
A simpler configuration consists of only the blade, the 
hub block and joint 1 all on a single load path. The 
control angles are then imposed as an angular motion 
command at the joint. The position and orientation of 
the joints are based on the UH-60A configuration. The 
properties of joint 1 — the stiffness and damping of the 
elastomeric bearing — are from the Army /NAS A master 
data base. The lag damper is assigned to this joint. The 
properties of joint 2 are unknown and left as zero. The 
pitch link stiffness is introduced either as a linear spring 
at the pitch link bottom joint (baseline configuration) or 
added to the pitch stiffness at the bearing (simple config- 
uration). Note that the elasticity of the pitch link adds 
to the net control system stiffness in the baseline case. 
For this reason the blade root pitch angles under trim 
conditions differ marginally for the two configurations. 

The baseline blade has an SC1095 profile (Fig 2) 
with a realistic internal structure. Each section (or seg- 
ment in 3D) consists of 37 hexahedral elements with 27 
nodes each (592 cross-sectional nodes in total) with 16 
spanwise segments in all. The grid is coarse and con- 
structed using a simple in-house grid generator. The 
cross-sectional grid is guided by Ref [15] but modified 
for second-order elements and then extruded in the span- 
wise direction accounting for quarter-chord axis, built- 
in twist, and tip sweep. A simpler blade has the same 
spanwise construction but a fake rectangular solid cross- 
section. The full-up model — baseline configuration and 
baseline blade — contains 17, 566 degrees of freedom. 

The baseline blade requires a detailed description of 
material properties. The material properties are unavail- 
able in the public domain so values are assigned such that 
they generate cross-sectional structural properties simi- 


Figure 3: Close-up of rotor hub. 

lar to the UH-60A. This requires a heterogeneous mate- 
rial assignment with different values for the box beam, 
skin, core, and leading edge elements. In general a 6 x 6 
material matrix can be assigned to each element; here 
isotropic properties are considered within each element. 
The same profile and cross-section is used at all span 
stations, thus, the chord length and structural properties 
remain uniform. The assignment of material properties 
to generate the target set of cross-sectional properties 
is carried out iteratively. The cross-sectional properties 
can be estimated from static deflections to prescribed 
forcing using the 3D analysis itself — a process similar 
to that of static testing of blades. 

For the extraction of sectional properties, and for 
the contraction of aerodynamic interfaces later, blade 
cross-sections must be defined. Figure 4 describes the 
definition (blade twist and some surface elements re- 
moved for illustration). Elements must supply nodes at 
the leading and trailing edges and these nodes have to 
be designated up-front as part of the model input. The 
analysis then defines the chord line and locates the 1/4— c 
(needed for CFD interface) and 3/4 — c (needed for clas- 
sical aerodynamics) points. Two additional nodes, one 
on the top surface and another at the bottom, are re- 
quired to define a nominal thickness line. It is used to 
simply define the cross-sectional plane. It need not be 
perpendicular to the chord line thus any two points suf- 
fice as long as they are on the top and bottom surfaces. 
The key requirement is that all element faces line up 
along the section. Thus even though the solver is natu- 
rally unstructured, sections can only be defined at span 
stations where the elements line up. The geometry and 
grid must be described accordingly for all span stations 
where sectional analysis is needed. 

If e c is an unit vector aligned along the chord line 


5 



Trailing 

edge 


Figure 4: Definition of cross-sectional planes. 



Figure 6: Calculation of elastic axis. 

(trailing edge to leading edge) and t an unit vector along 
the thickness line (bottom surface to top surface) then 
the unit vectors along the axial and normal directions 
are defined as the cross products 
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As the blade deforms the orientation of e a , e c , e„ can be 
used to extract a set of beam-like cross sectional rota- 
tions. Elastic deformations can be found by subtract- 
ing the baseline geometry consisting of the built-in twist 
(swept tip consists of sheared sections for this rotor). 
The deflection of the 1/4 — c point can be used as the 
cross-sectional linear deflection. 

The sectional property calculation is performed as 
follows. Note that this step is not needed for 3D analysis 
but only used as a means for backing out the unknown 
3D material properties. There is no unique inverse solu- 
tion so an approximate representation is sought by trial 
and error. The inertial properties and flexural constants 
are obtained directly from geometry and material de- 



Figure 5: Calculation of torsional rigidity GJ. 
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Figure 7: Cross-section of the main blade. 

scriptions. The torsional stiffness GJ and the shear cen- 
ter require static solves. The steps are as follows: 

1. Obtain sectional mass, center of gravity (C.G.), area 
and center of area (C.A.), and axial stiffness EA. 

2. Obtain Principal Axes, bending stiffnesses flap and 
lag El about C.A. along Principal Axes, moments of 
inertia about C.G. along Principal Axes, and C.G. 
and C.A. offsets from 1/4 — c along Principal Axes. 

3. Introduce a pure couple at the blade tip; calculate 
spanwise twist gradient < j>' and reaction torque r by 
stress integration. Calculate GJ = t/4>' . 

4. Introduce a series of forcing at the tip in vertical and 
chord wise directions to identify the shear center. 

Steps 1 and 2 are straight forward. For steps 3 and 4 
an uniform blade is constructed (zero built-in twist and 
sweep) with the cross section. In step 3 a pure couple is 
introduced at the blade tip using two equal and opposite 
nodal forces along any direction. The value of GJ calcu- 
lated at sections away from the ends remain uniform as 
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shown in Fig 5 and considered to be the cross-sectional 
value. The calculations deviate near the root due to edge 
effects and the tip due to proximity to applied forcing as 
expected. For the rectangular cross-section the classical 
St. Venant constant is recovered. The calculation of the 
shear center is show in Fig 6. Elastic twist from a set 
of vertical and chordwise forcing are interpolated to find 
the zero cross over point and identify the location of the 
shear center. For no bending-twist coupling the shear 
center is a cross sectional property and remains uniform 
with span. 

The blade cross-section with the calculated offsets 
are shown in Fig 7. Note that the blade is placed with 
pitch axis ahead of the rotation axis by the torque offset. 
The C.G. is placed close to and ahead of the pitch axis 
at 1/4 — c. Then stability is guaranteed in the thin airfoil 
sense (aerodynamic center at 1/4 — c). But the E.A. falls 
significantly behind the C.G. (by 13% c) which results in 
a strong flap-torsion coupling. The effect of this coupling 
is seen later in airloads. The C.G. location relative pitch 
axis is similar to that of the UH-60A but its offset relative 
E.A. is a substantial deviation. 

3D CFD/CSD FORMULATION 

The 3D CFD/CSD formulation includes the inter- 
face (spatial coupling) and solution procedure (temporal 
coupling) . 

The interface deals with deformations sent to CFD 
and surface forcing sent to CSD. The deformations can 
be described in one of two ways: as a 2D beam-like field 
or as a 3D field. The surface forcing can be described 
in one of three ways: as 2D segmental airloads, as 3D 
patch forces, or as 3D surface pressures and shear. It 
is assumed that the discretized surface will always dif- 
fer between CFD and CSD and hence corrections will be 
required for conservation (same virtual work) and con- 
sistency (same integrated forcing). The difference is not 
merely from mesh mis-match but unavoidable geometry 
mis-match — aerodynamic design will define the sur- 
face first while structural design will populate the sur- 
face last. Thus a method is desired that can be applied 
even when the internal structure is at an early design 
stage and does not necessarily carry on to the wetted 
surface. 

Deformation Interface 

A 2D beam-like description of the deformation field 
allows the airfoil shapes in the fluid domain to be kept 
intact. While this description is the natural state of the 
art in beams, in 3D, they must be extracted from the 
deflection field. The extraction is based on the deformed 
chord and thickness lines that describe the deformed ori- 
entation of the cross sectional plane. The cross sectional 
planes are defined as shown earlier in Fig 4. The extrac- 
tion excludes any cross-sectional flexibility — beyond 


what is captured by the four nodes that define the lines 
— and appears to defeat the purpose of a 3D interface 
but is indeed the only means by which structures that are 
internally incomplete can be analyzed. Thus integrated 
3D stress/strains can still be obtained at early stages 
of structural design when only the main load-bearing 
pieces may be in place. The extraction procedure is 
straight-forward: the direction cosines of the deformed 
cross-section: e Q ,e ra ,e c , are used to extract three Euler 
angles for any order of rotation. Standard limitations 
near ±7r/2 apply but seldom encountered in rotor prob- 
lems. The order used here is lead-lag, flap, and pitch, 
consistent with mesh deformation. 

A 3D description of the deformation field encounters 
the problem of geometry and mesh mis-match. Ideally, 
there is a common 3D CAD geometry that populates 
both sides of the interface — CFD and CSD — so that 
a point in either domain can be associated with a cor- 
responding point on the other via geometry. A common 
geometry description is beyond state of the art so the 
CFD mesh is presently considered to be the exact ge- 
ometry. The problem then reduces to associating each 
CFD mesh point with its counterpart on the CSD sur- 
face. This parameterization is different in 3D and is 
the basis for constructing an exact interface. Here ex- 
actness is defined as errors introduced due to interface 
mis-match not exceeding those that are introduced by 
the discretization of the individual solvers themselves. 
The parameterization is described under Level-11 Force 
Interface later. 

Level-I Force Interface 

The level-I force interface is an extension of the 
lower order aerodynamic interface with 2D airloads ad- 
mitted from CFD. Guaranteeing conservation is not pos- 
sible but consistency is achieved by the use of segmen- 
tal airloads (integrated over span-wise segments, see 
Refs [14, 17] for methodology and validation). These 
airloads must now be distributed over the surface nodes 
of the segmental 3D elements. The nature of 2D air- 
loads implies the 3D finite element grid must ensure sets 
of surface elements along spanwise strips. Elements can 
be unstructured within the strip but set boundaries must 
line up along the chordwise direction. Then all surface 
nodes within a segment are identified as aerodynamic 
interface nodes. 

In the present mesh the 3D elements are naturally 
aligned in the spanwise direction therefore segments are 
easy to define. Figure 8 shows two spanwise element 
rows (blade twist and some surface elements removed 
for illustration). Because the elements carry internal 
nodes there are five spanwise nodal lines. Assume the 
two ends represent the root and the tip. Then the aero- 
dynamic segments span half-way across either side of 
each nodal line. So there are as many segments as there 
are nodal lines. Each segment receives dimensional air- 
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Figure 8 : Level-I force interface. 

loads from CFD with pitching moment about the 1/4 — c. 
These are then distributed over the segmental nodes in 
an assumed manner — linearly for lift (along chord, in- 
creasing toward leading edge) and quadratically for drag 
(along thickness) — as functions of three constants to be 
solved at each azimuth and radial station from the CFD 
airloads. If y and z are the chordwise (toward leading 
edge) and thickness wise (to upper surface) coordinates 
and Fz, Fy and M are the vertical, in-plane and 1/4 — c 
pitching moments then the nodal forces fz and fy in 
the vertical and in-plane directions are represented as 

fz = a 0 + a 1 (y-y 25 ) 
fy = ( 3 o + ai (z - z 25 ) 2 

with the three constants 00,01 an d (3o solved using the 
three airloads 

Fz = ^ fz 

Fy = J2fy (3) 

M 2 5 = ^2 ifz(y ~ 2/25) — fy{z ~ 225)] 

The subscript 25 denotes 1/4 — c quantities. The sum- 
mation is over all n segmental nodes. The axial force 
Fx, although negligible compared to the inertial force, 
is simply distributed at each node as fx = Fx /n. 

The internal stresses from the level-I interface are 
only as good as the assumed representation and there- 
fore ad hoc. It is a convenient interface for initial veri- 
fication because it leaves the Delta Coupling procedure 
intact. This interface is used later for a fully coupled 
solution. But a 3D CSD model opens opportunity for a 
detailed interface with which exact internal stresses can 
be calculated. This is described below. 

Level-II Force Interface 

The level-II force interface is an exact 3D patch force 
interface. Consistency is guaranteed but conservation 
achieved only at the limit of fluid and structural mesh 


refinements. In this interface each CFD surface mesh 
point (or the center of each surface triad or quad — 
depending on how the surface stresses are integrated in 
the fluid domain to generate a patch force) is first asso- 
ciated uniquely with a CSD surface coordinate. Ideally, 
as noted earlier, a point on the exact geometry should be 
associated not a point on the CFD mesh, but in absence 
of a geometry module that populates both sides of the 
interface, the CFD mesh is assumed to be the exact ge- 
ometry. The method is generic to mis-matched meshes. 
Given a CFD surface point — regardless of whether it 
lies on the CSD surface or not — an unique set of CSD 
surface parameters are calculated. This parameteriza- 
tion is based only on the natural coordinates and shape 
functions of the elements. 

Within an element (including its surface) the geo- 
metric coordinates are related to the curvilinear natural 
coordinates as 

x = x(£, Hi C) 

y = y(t,v,0 ( 4 ) 

« = z(f,»7, C) 

where —1 < £, 77 , / < 1. Because the relationship is non- 
linear, given a target point (x , y , z)t, its corresponding 
natural coordinates (£ 77 f) must be found iteratively. 
The iterations can be started from any of the element 
nodes — where both (£ 77 </)o and ( x y z)o are known — 
and continued as per Newton-Raphson to convergence. 

k = 0,1,2 ,... 
x v 

Vi Vv vc 

Z n Z( 

( 5 ) 

The left hand side gradient matrix (where x £ = dx/df, 
and so on) is the transpose of the Jacobian of the co- 
ordinate transformation. As long as the target point 
lies within or on the element and the element is well- 
behaved (i.e. has a well-conditioned Jacobian) the iter- 
ations will converge to an unique solution. Thus if the 
CFD point is on the surface (matched interface) its nat- 
ural coordinates can be readily obtained using this basic 
procedure. A generic procedure is one that will obtain 
a corresponding set of natural coordinates for all points 
including those that are out of surface but collapse to 
the basic procedure for points on surface. Two different 
generic procedures are proposed as follows. 

Method 1: Gradient extrapolation. Let £ = 1 be the 
surface (see Fig 9). An out of surface point is associ- 
ated with a surface point (£, 77 ) if the surface coordinate 
£ leaves the surface in such a direction that it passes 
through the point. Thus the method can be termed a 
gradient extrapolation method as the gradient of the sur- 
face coordinate (£) is extrapolated. Note that all out of 
surface points that lie along the direction of extrapo- 
lation will be associated with the same — but unique 



— CSD point. To identify this point the same itera- 
tions as Eq 5 are used but the gradient matrix is always 
evaluated at the surface (£ — 1 in this example). The 
iterations will drive the absolute value of £ to greater 
than 1, confirming the point lies outside the surface, but 
will converge in £ and rj. If not, or if the absolute values 
of either of these converged coordinates are greater than 
1, then the point is rejected by the £ = 1 surface. All 
surfaces defined as fluid interfaces are checked one by 
one and if the point is rejected by all it is rejected by the 
element. If rejected by all elements it is eliminated from 
the interface and tagged instead as a candidate for er- 
ror calculation. Figure 9 shows an example of successful 
extraction of surface parameters using this method. 


Method 2: Orthogonal projection. Let £ = 1 be the 
surface (see Fig 9). An out of surface point is associated 
with the nearest surface coordinate (£, rf) (note that this 
is not the nearest neighbor method where the association 
is to the nearest node). This means the normal to the 
surface at that coordinate passes through the point. The 
surface coordinates are obtained by 


km 0,1,2,... 


x v 

V( Vr, 
Z£ Zr, 




A Moore-Penrose inverse of the left hand side matrix 
(A + = (A 2 A)~ 1 A r ) solves the above equation in the 
least square sense and finds the nearest surface coordi- 
nate. If the absolute values of either of the converged 
coordinates are greater than 1 then the point does not 
belong to the £ = 1 surface. As in method 1 all surfaces 
of all elements are checked one by one for inclusion or 
elimination from the interface. Figure 9 shows an exam- 
ple of successful extraction of surface parameters using 
this method. Figure 10 shows an example of rejection 
where the extracted surface parameters lie outside the 
element. 

The surface coordinates produced by the two meth- 
ods are different (Fig. 9). In general the farther the point 
from the surface the greater the difference. Only when 
the element boundaries are orthogonal the two methods 
converge. Because then the coordinate leaving the sur- 
face (used in method 1) is the surface normal (used in 
method 2). In general either method can be used with 
the difference compensated for by the interface error cal- 
culation (see below). For each CFD surface node the 
corresponding CSD element and its natural coordinates 
£, r ], £ then complete the surface parameterization. The 
displacement of each CFD node and the virtual work 
contribution of the patch force occurring at that node 
are now precisely defined. 


Force interface error. Even though the interface is 
exact — the errors now stem entirely from geometry de- 
scription of the individual solvers and not how variables 


CFD 



Figure 9: An example of an out of surface CFD 
point for which CSD surface parameters could 
be extracted successfully. 


CFD 

surface 

Iterative search for point 

CSD surface point 





Figure 10: An example of an out of surface 
CFD point for which CSD surface parameters 
could not be extracted successfully. 


are exchanged between them — it can lead to significant 
errors or can even break down (large number of points 
eliminated from interface) when the underlying struc- 
tural model is incomplete. It is important that these 
errors be accounted for so that meaningful stress calcula- 
tions are still possible at early stages of structural design. 
For this purpose the error is defined in terms of segmen- 
tal airloads. For all CFD nodes that are eliminated from 
the interface (i.e. , not claimed by any CSD element) the 
error is simply its contribution to the segmental airload. 
For all others, the error is only in form of a segmental 
pitching moment, calculated based on the distance be- 
tween the CFD coordinate and its corresponding CSD 
coordinate. The segmental airloads (errors) are then re- 
introduced into the structure through a level-I interface. 
For no mis-match the error vanishes. For a 100% mis- 
match level-II collapses to level-I — since every point is 
then rejected and re-introduced through level-I. 
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Figure 11: CFD and CSD surface meshes for an UH-60A rotor. 


UH-60A example. The UH-60A near body CFD 
mesh is superposed on the UF[-60A-like CSD surface el- 
ements in Fig 11. As expected there are significant mis- 
matches in geometry and meshing. The surface parame- 
ters are now extracted using the two methods described 
earlier. Because the structural model has a coarse sur- 
face mesh a simple search algorithm suffices — each 
point is searched over all surface elements until claimed. 
The search process requires an additional check. For ex- 
ample if C = 1 is the assigned surface then a CFD point 
is parameterized for both surfaces £ = ±1 and is ad- 
mitted only if it is closer to the assigned surface. This 
prevents a CSD lower surface element claiming a point 
from the CFD upper surface. Pathological cases are still 
possible however for dramatic mis-match. For example 
if the CFD mesh is entirely below the CSD mesh then 
all points are claimed by elements at the lower surface 
and none by elements at the top. A common geometry 
description is needed to prevent such pathological cases. 

Figures 12 and 13 show the eliminated points from 
methods 1 and 2 respectively. Most of these elimina- 
tions are easily explained — the structural model has no 
tab and there is a discrepancy in chord length on the 
swept tip portion — but method 2 eliminates a signifi- 
cant number of points near the leading edge. One such 
point is examined closely in Figs 14 and 15. Figure 14 
shows that the point is eliminated because none of the 
elements provide a surface normal that passes through 
that point. The reason is that 3D elements are only 
C° (only displacements are continuous not slopes) thus 
there is a discontinuity in surface normals between ele- 
ments. Points that lie in that gap are hidden from the 
interface. This is avoided by the gradient extrapolation 
method as shown in Figure 15. The gradient of the out of 


surface coordinate is always continuous across elements 
as long as they share a common face. Thus if the struc- 
tural mesh ensures surface elements with common faces 
then the gradient extrapolation method appears to be 
the preferred method. 



Figure 12: CFD surface points (in red) elimi- 
nated from interface by the gradient extrapo- 
lation method (method 1). 


Level-Ill Force Interface 

The level-III interface is a 3D surface pressure and 
shear force interface. The expression for virtual work 


10 


Element 1 



Assigned points 23712(95.8%) 
Unassigned points 1049 (4.2%) 


Total points 24761 


Figure 13: CFD surface points (in red) elimi- 
nated from interface by the orthogonal pro- 
jection method (method 2). 



Figure 14: A leading edge CFD surface 

point eliminated by the orthogonal projection 
method. 


integrates terms involving both fluid stresses and struc- 
tural surface normals. Because conservation requires 
all interpolation to be carried out using corresponding 
shape functions in each domain [16], either fluid stresses 
must be supplied at all structural Gauss points (if in- 
tegrating in the CSD domain) or structural shape func- 
tions supplied at fluid mesh points (if integrating in the 
CFD domain). Then conservation is guaranteed but con- 
sistency is still achieved only at the limit of fluid and 
structural mesh refinement. For different mesh resolu- 
tions the integrated forces will remain different in CFD 
and CSD domains. This is undesirable for rotor prob- 
lems where accuracy of the trim state is essential. Thus 
the level-II interface is considered the most desired. The 
level-I interface is an approximation of level-II and is 
used for integrated analysis as the first step. 


CFD point 



Figure 15: A leading edge CFD surface 

point admitted by the gradient extrapolation 
method. 


CFD/CSD solution procedure 

The level-I force interface is used. Thus the con- 
ventional Delta Coupling procedure of Johnson (Loose 
Coupling in rotorcraft terminology) can be retained with 
segmental airloads (dimensional) used as the delta vari- 
ables. (Ref [18] is the original invention, the current im- 
plementation follows Ref [19]). The level-II force inter- 
face requires an advanced formulation with the surface 
element generalized loading used as the delta variables. 
This method was proposed and validated in Ref [17] for 
beams but not yet implemented in the current solver. 
The 2D deformation interface is used. Thus the conven- 
tional beam-like mesh deformation procedures can be 
retained in the CFD solver. The 3D deformation inter- 
face requires a mesh deformation method that is beyond 
scope of the current solver. Thus an integrated solution 
is obtained using the simplest of the interfaces. 

As per the Delta Coupling procedure periodic air- 
loads and deformations are exchanged. This accommo- 
dates unequal time steps in the fluid and structural do- 
mains naturally — 0.25° and 3° respectively — with 
airloads and deformations digitally filtered (Fourier in- 
terpolation, 12 harmonics used here) during exchange. 
Although airloads are available at the CSD time steps, 
deformations are not available at the CFD time steps, 
hence an interpolation is always needed. But because the 
CSD solution is obtained by time integration this inter- 
polation cannot be made consistent with the solver while 
at the same time provide smooth grid motion across CFD 
time steps. Thus time accuracy cannot be ensured — 
not even at 3° intervals. The alternative — to march 
CFD and CSD together with intermittent sub-iterations 
— though straight forward is practically unacceptable in 
rotorcraft as the structural dynamics takes 30 — 40 revo- 
lutions to attain the trimmed periodic solution whereas 
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Figure 16: Flap moment at root bearing. 



Figure 18: Pitch moment at root bearing. 


the fluid solution settles rapidly to periodicity within 
1 — 2 revolutions. And because the fluid solution set- 
tles rapidly, advancing the solver by only a quarter of 
a revolution can allow nominally periodic airloads to be 
constructed using all four blades. The coupling process 
uses half a revolution for the first two CFD iterations 
and quarter revolutions there after. 

As per the Delta Coupling procedure the simple 
internal aerodynamic model is retained in CSD as a 
preconditioner for convergence. An isolated rotor trim 
model is used with thrust and two hub moments targeted 
using collective and cyclic commands imposed as pitch 
link displacements. The trim Jacobian is calculated once 
using the internal aerodynamic model and requires no 
update for all coupling iterations for the flight condition 
studied. 

3D CFD/CSD ANALYSIS RESULTS 

The UH-60A flight test counter 8534 (advance ra- 
tio fj, = 0.368, nondimensional thrust CV/cr = 0.084 and 
shaft tilt as = —7.31° (forward)) is considered. It is 
a high speed (158 kts) high vibration level flight where 



Figure 17: Lead-lag moment at root bearing. 



Figure 19: Joint reactions at root bearing. 


the dominant airloads are the 3D unsteady tip transonic 
pitching moments and hence CFD is essential. Even 
though the blade structural properties are not identical 
to UH-60A the aerodynamic characteristics (geometry 
and grid) are the same and therefore comparing predic- 
tions with UH-60A test data provides a reasonable basis 
for qualitative validation. 

Dynamic Response 

The time step required for stability and convergence 
of the dynamic response is observed to be governed by 
the nonlinearities introduced by the joint not the size 
of the finite element mesh. Thus the simple blade and 
the baseline blade both require the same time step for 
the same hub configuration. It is also observed that the 
dominant source of nonlinearity is the elastomeric bear- 
ing where the trim commands are imposed. Hence the 
simple hub and the baseline hub also have the same time 
step requirement — as the bearing is common to both. 
The strong nonlinearities require the calculation of dy- 
namic stiffness at every time step. Calculating it once 
about a trim solution or updating it every few time steps 


12 


diverges the solution quickly. Under these considera- 
tions both time integration methods (the Generalized-a 
was executed at its simplest Newmark form) converges 
to the same results. Sub-iterations are possible at ev- 
ery time step but needed only for relatively large steps. 
Typically 3 — 5 sub-iterations were required to drive the 
residual down to 10 -5 at each time step. Sub-iterations 
are also required for the consistency of inertial loads with 
the response solution. For the purposes of CFD coupling 
smaller time steps were preferred over larger steps with 
sub-iterations so as to minimize errors from deformation 
interpolation. Thus instead of A tp = 9° with 3 — 4 sub- 
iterations, Ail’ — 3° with no sub-iterations was preferred. 
The results shown are all obtained using Ai/j = 3° with 
no sub-iterations. 

The blade dynamic response is verified by studying 
the inertial loads and joint reactions. Simple aerody- 
namics (without CFD) is adequate for this purpose. To 
verify root reactions the simpler configuration is used 
where the pitch link stiffness is lumped at the bearing 
and the control angles are imposed as command sig- 
nals at the bearing. Figures 16, 17 and 18 show the 
integrated aerodynamic, inertial and total (aerodynamic 
plus inertial) moments at the elastomeric bearing for flap 
(positive down), lead-lag (positive forward) and pitch 
(nose up) motions. The equal and opposite nature of 
the aerodynamic and inertial flap moments verifies the 
fundamental dynamics of the blade. It is not identi- 
cally zero because of a small amount of stiffness and 
damping at the bearing. The damping is significantly 
larger in lag (due to the lag damper) which is reflected 
in the more substantial lead lag moment. The stiffness 
is significantly larger in pitch (due to pitch link stiffness) 
so that the root reaction follows the aerodynamic forc- 
ing closely. The inertial loading calculation, and hence 
the velocity and acceleration calculations, are verified by 
comparing these integrated moments obtained by force 
summation with reaction forces obtained directly from 
joint response at the bearing. This comparison is af- 
fected by the solver time step as mentioned earlier but 
agrees closely for A^ = 3° as shown in Fig 19. The joint 
reactions are left in the joint frame to illustrate the effect 
of order of rotations. The force summation results are 
along the undeformed blade frame. Because the order of 
joint rotation is lead-lag first then flap and then torsion 
only the lead-lag comparison is exact, the flap and tor- 
sion moments have to be transferred to the undeformed 
frame for exact verification. Henceforth all results shown 
use the force summation method. 

CFD/CSD Airloads 

The convergence of the CFD/CSD analysis proceeds 
as shown by the mean normal force distribution from 
CFD in Fig 20. Iteration-0 shows CFD airloads calcu- 
lated using deformations from the baseline trim solu- 
tion (with internal aerodynamics). By Iteration-7 the 


airloads and dynamic response have nearly converged. 
Henceforth all airloads are shown from both the last two 
iterations 6 and 7. It is known that the flight test normal 
force shows a greater thrust at this condition than what 
can be achieved by the trim target. 



Figure 20: Steady normal force convergence from 
3D CFD/CSD using level-1 interface; predictions 
compared with measured UH-60A airloads for 
qualitative comparison. 

The detailed sectional airloads at three span wise 
stations are shown in Fig 21. Clearly the normal forces 
deviate significantly from the UH-60A values particu- 
larly at the outboard stations (inboard of 67.5%R predic- 
tions are good and hence not shown). The discrepancy is 
partly because of the less severe transonic pitching mo- 
ment drop at the first quadrant close to the tip (96.5%R) 
due to a different blade dynamics than the UH-60A but 
mainly for reasons revealed in Fig 22. Figure 22 shows 
the harmonics of the normal force across the blade span. 
Forcing that are passed to CSD is plotted (divided by 
segment lengths) so the spanwise resolution is the same 
as the number of aerodynamic segments in the model (33 
segments). The 1/rev phase with its distinct 180° shift 
over the span verifies the accuracy of the trim solution. 
All the vibratory harmonics (3—5/rev) are well predicted 
with at least the same accuracy achieved by current UH- 
60A beam models coupled to CFD. The main discrep- 
ancy is from the 2/rev airloads. These are determined by 
a large 1/rev torsion at this high speed condition. The 
1 / rev torsion in the current analysis deviate significantly 
from the UH-60A due to the very high (13% c) C.G. off- 
set from the E.A. (Fig 7). The 1/rev torsion magnitude 
is as expected (about negative 4° peak to peak) but the 
phase is contaminated by 1/rev flap. Thus a proper val- 
idation of a 3D analysis requires a more careful assign- 
ment of material properties so that at least the major 
offsets (E.A. and C.G.) are precisely placed. 
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Figure 21: Airloads from 3D CFD/CSD using level-1 interface; predictions from last two iterations 
compared with measured UH-60A airloads for qualitative comparison. 
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Figure 22: Normal force harmonics from 3D 

CFD/CSD using level-1 interface; predictions 
compared with UH-60A airloads for qualitative 
comparison; symbols are flight test measure- 
ments, lines are predictions from last two iter- 
ations. 


CFD/CSD Dynamic Stresses 

A 3D analysis enables the direct calculation of dy- 
namic stresses and strains on all 3D flexible parts. Re- 
sults from the blade are shown here. There is no data 
readily available for validation but because blade loads 
(bending moments) are measured by calibrating surface 
strains, validation is possible in principle. 

The predicted tensile stresses (Cauchy stresses r n 
where 1 is nominally the radially outward direction) oc- 
curring at 77% R are shown in Fig 23 for every 45° az- 
imuth. The airfoil section (expanded in scale) is over- 
layed for reference. The mean stress level is set by the 
centrifugal loading; oscillatory stresses arise from bend- 
ing. As expected it is the main spar — the box beam 
at the center — that has the maximum loading and en- 
counters the maximum variation. At 90° azimuth the 
maximum stresses occur at the flange-like top portion 
of the spar whereas at 270° they occur at the bottom. 


In general the top portion is loaded more heavily on the 
advancing side and the bottom portion on the retreating 
side. The loading at 0° and 225° are similar on the top 
and bottom implying that the bending curvatures are 
possibly close to zero at these azimuths. The web-like 
vertical portions of the spar do not show any significant 
oscillatory loading. 

The loading pattern is studied in greater detail at 
eight locations on the cross-section (Fig 24) — two on 
the trailing-edge (one each near the top and bottom sur- 
faces), four inside the spar (one on each arm) and two 
near the leading-edge (on each near the top and bottom 
surfaces). Figure 24(a) shows these sensor locations. In- 
stead of stresses micro-strains (fi) are plotted versus az- 
imuth (Figs 24(b)- 24(d)) since maximum safety of flight 
limits are often characterized by strains. These reflect 
not only the nature of the loading but also indicate the 
extent of cross-sectional flexibility. The strains are sig- 
nificant all over the cross-section. At the trailing-edge 
(Fig 24(b)) — where the sensors are close together — 
the strains are similar on the top and bottom surfaces as 
expected. But on the spar they exhibit a pattern similar 
to flap bending moments (Fig 24(c)). The typical rise in 
flap bending moment in the third quadrant is reflected 
as a reduction in strain on the top (increasing compres- 
sion) and an increase in strain at the bottom (increasing 
tension). The pattern is almost symmetric except for an 
additional effect from chordwise bending. The peak-to- 
peak oscillatory strains are around 800/1 on the top and 
bottom portions of the spar. The web-like portions on 
the back and front encounter relatively benign strains. 
The strains here reflect the chord wise bending moments. 
The chord wise stiffness is significantly higher than the 
flap wise stiffness and this is reflected in the lower strains. 
The strains at the leading-edge (Fig 24(d)) also exhibit 
similar characteristics as flap bending but the location of 
the top sensor is such that its behavior is closer to that 
of the front web sensor rather than the top surface. 

The transverse shear stresses in the vertical direc- 
tion (Cauchy stresses T13 where 1 and 3 are nominally to- 
ward the tip and the top surface respectively) are shown 
in Fig 25. In general the shear stresses are an order of 
magnitude lower than the tensile stresses. Unlike ten- 
sile stress here the web-like vertical portions carry the 
maximum loading. Between 45°-180° the front web is 
more heavily loaded whereas from 180°-45° the loading 
is transferred to the back. But in general the rear web is 
loaded more heavily reflecting the generally nose down 
characteristics of the torsion loading — except in be- 
tween 90°-135° where it possibly reflects the local wake 
induced loading in the advancing side. The rear web has 
generally greater stress levels than the front web per- 
haps affected by the relative positions of shear center 
and aerodynamic center. The transverse shear stresses 
in the in-plane direction (Cauchy stresses T12 where 2 is 
nominally toward the leading-edge) are shown in Fig 26. 
These are similar in magnitude to those in the vertical 
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(e) Azimuth -0 = 180° 


(f) Azimuth ?/> = 225° 



Figure 23: Dynamic tensile stresses (rn) predicted using 3D CFD/CSD analysis; 77% R station; 
colormap range: 10 6 — 18 x 10' N/m 2 (Pascal). 
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(a) Strain sensor locations; two at the trailing-edge, four on the main 
spar and two near the leading-edge. 
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(b) Strains at the trailing-edge. 
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(d) Strains at the leading-edge 


Figure 24: Azimuthal variation of tensile strains (en) predicted using 3D CFD/CSD analysis; 77% R 
station. 


direction but they occur at different locations. Like the 
tensile stresses, here the flange-like portions of the spar 
exhibit the maximum loading. But in terms of gross 
loading they exhibit the same general pattern as the 
vertical shears — the rear web is loaded more heavily 
reflecting the generally nose down characteristics of the 
torsion loading except in between 90°-135° azimuths. 

The generic nature of the blade and the approximate 
nature of the level-I interface make specific conclusions 
on the details of stress/strain variations premature for 
this rotor but the results above appear consistent with 
typical patterns of loading and reveal the type of detailed 
stress/strain analysis possible using 3D structures. 

CONCLUSIONS 

A full three-dimensional finite element-multibody 
structural dynamics solver was coupled to a three- 
dimensional Reynolds-averaged Navier-Stokes solver for 
the prediction of integrated stresses and strains on a ro- 
tor blade in forward flight. All major pieces of the 3D 
workflow were implemented and analyzed. The struc- 


tural model was an idealized representation of a fully 
articulated UH-60A-like rotor. The aerodynamic model 
was identical to the UH-60A rotor. The primary focus 
of this work was on the 3D CFD/CSD coupled solution 
procedure. Two different types of 3D interfaces were de- 
veloped and the simpler of the two implemented. The 
coupled airloads were compared with flight test data. 
The coupled 3D blade stresses and strains were stud- 
ied for generic patterns typical of rotors in high speed 
forward flight. Based on this work the following key 
conclusions are drawn. 

1. An integrated aeromechanical analysis with a high- 
fidelity 3D representation of the structure can in- 
deed be carried out for helicopter rotors in forward 
flight. Dynamic stresses and strains resulting from 
aerodynamic loading can then be calculated with- 
out any reduced order approximation. Advanced 
hub types and blade shapes can then be modeled 
without any reduced order assumption. 

2. The convergence of the (implicit) CSD solver re- 
quires dynamic stiffness to be updated at every time 
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(g) Azimuth 0 = 270° 


(h) Azimuth 0 = 315° 


Figure 25: Dynamic shear stresses in vertical direction (T 13 ) predicted by 3D CFD/CSD analysis; 77% 
R station; colormap range: —5 x 10 6 — 20 x 10 6 N/m 2 (Pascal). 
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(a) Azimuth -0 = 0° 



(c) Azimuth 0 = 90° 




(b) Azimuth 0 = 45° 



(d) Azimuth 0 = 135° 



(e) Azimuth 0 = 180° 



(f) Azimuth 0 = 225° 
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(g) Azimuth 0 = 270° 


(h) Azimuth 0 = 315° 


Figure 26: Dynamic shear stresses in in-plane direction (ti 2 ) predicted by 3D CFD/CSD analysis; 77% 
R station; colormap range: —5 x 10 6 — 20 x 10 6 N/m 2 (Pascal). 
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step. This is due to the strong non-linearities intro- 
duced by joints unified within 3D elements. For rel- 
atively large time steps sub-iterations are required 
— both for the accuracy of the response solution 
as well as for consistency of inertial loads with the 
response solution. 

3. The time step is not dictated by the mesh resolu- 
tion of the flexible 3D parts but rather by the non- 
linearities introduced by the joints. Thus a fine- 
mesh problem could be solved with the same stabil- 
ity and convergence behavior using the same time 
step as long as the trim conditions (joint motions) 
remained similar. 

4. it is not essential that the finite element internal 
structure be populated in every detail for a nominal 
application of the analysis. The 3D interface can 
be made generic enough to accommodate an incom- 
plete internal structure provided the primary load 
bearing pieces are in place. But a careful descrip- 
tion of the material properties of these pieces are 
required so that the major blade offsets (E.A. and 
C.G.) are correctly located. 

In the end we note that even though the CSD solver 
is parallel and scalable it was executed on a single pro- 
cessor in this study. The structural geometry, mesh, 
and material descriptions were representative of a re- 
alistic rotor but inexact and coarse. The exact patch 
force based level-II interface was developed but only the 
approximate level-I interface was implemented for inte- 
grated solution. A time integration procedure is an inef- 
ficient solution process to reach a periodic trim solution 
in a low-damped near resonance system like rotorcraft. 
These and other topics remain the subjects of future 
work. 
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